!----------------------------------------------------------------------------
! Program: 2D_Potential_FMM - A fast multipole boundary element
! Method (BEM) code for analyzing large-scale, general 2D
! potential problems (governed by Laplace equation)
! using constant elements.
!
! Developers: Dr. Naoshi Nishimura at Kyoto University, Japan;
! Dr. Yijun Liu at the University of Cincinnati, Cincinnati, OH.
!
! Version: V.1.20.
! Released: October 1, 2008.
!
! Copyright(c)2004-2008 By Kyoto University and University of Cincinnati.
! This code is intended for educational use only. No part of
! the code can be used for any commercial applications/
! distributions without prior written permissions of the
! original developers.
!
!----------------------------------------------------------------------------

program fmm
   use kind_parameter
! Purpose:
! To illustrate some of the basic features of a Fortran program.

   integer:: maxl, levmx, nexp, ntylr, maxia,&
      ncellmx, nleafmx,mxl,nwksz,&
      n,nfield
   integer :: io_in_fmm, io_in_dat, io_out, io_phi_bnd, io_xy, io_phi_dmn
   real(dp)::tolerance
   character(len=80) Prob_Title
   integer, allocatable:: ia(:)
   complex(dp), allocatable :: am(:)


   call CPU_Time(time0)
   open (io_in_fmm, file='input.fmm', status='old')
   open (io_in_dat, file='input.dat', status='old')
   open (io_out, file='output.dat', status='unknown')
   open (io_phi_bnd, file='phi_boundary.plt', status='unknown')
   open (io_xy, file='xy.plt', status='unknown')
   open (io_phi_dmn, file='phi_domain.plt', status='unknown')

! Input the parameters
   read(io_in_fmm,*)   maxl, levmx, nexp, ntylr, tolerance
   read(io_in_fmm,*)   maxia, ncellmx, nleafmx, mxl, nwksz
   read(io_in_dat,'(a80)')   Prob_Title
   read(io_in_dat,*)   n, nfield
   write(io_out,'(a80)') Prob_Title
   write(*,'(a80)') Prob_Title

! Estimate the maximum numbers of the cells and leaves,
! and size of the preconditioning matrix, etc.
   if(ncellmx.le.0) ncellmx = max(4*n/maxl,100)
   if(nleafmx.le.0) nleafmx = max(ncellmx/2,100)
   if(nwksz .le. 0)then
      nwksz = maxl*maxl*nleafmx
      ligw = 20
      lrgw = 1+n*(mxl+6)+mxl*(mxl+3)
      iwksz = n+3*nleafmx+1
      allocate (ia(maxia))

! Load the addresses (pointers) associated with the locations of the
! variables to be stored in the large array "am"
      call lpointer(lp, ln, maxia, ia, n, nexp, ntylr, ncellmx,&
         levmx, ligw, lrgw, nwksz, iwksz, nfield,&
         l_n, l_x, l_y, l_node, l_dnorm,&
         l_bc, l_a, l_b, l_xmax,&
         l_xmin, l_ymax, l_ymin, l_ielem, l_itree,&
         l_level, l_loct, l_numt, l_ifath, l_lowlev,&
         l_maxl, l_levmx, l_nexp, l_ntylr, l_tolerance,&
         l_ncellmx, l_nleafmx, l_mxl, l_u, l_ax,&
         l_sb, l_sx, l_ligw, l_lrgw, l_igwk,&
         l_rgwk, l_nwksz, l_iwksz, l_rwork, l_iwork,&
         l_xfield, l_nfield, l_f)

! Estimate the memory usage

      maxa = lp
      write(io_out,100) maxa*16/1.D6
      write(*,100) maxa*16/1.D6
100   format(/'Memory size of the large block am =', f12.1,' Mb'/)

! Allocate the large block 'am'
      allocate (am(maxa))

! Assign the parameters to the array am()
      call assigni(n, am(l_n))
      call assigni(maxl, am(l_maxl))
      call assigni(levmx, am(l_levmx ))
      call assigni(nexp, am(l_nexp))
      call assigni(ntylr, am(l_ntylr))
      call assignd(tolerance, am(l_tolerance))
      call assigni(ncellmx, am(l_ncellmx))
      call assigni(nleafmx, am(l_nleafmx))
      call assigni(mxl, am(l_mxl))
      call assigni(ligw, am(l_ligw))
      call assigni(lrgw, am(l_lrgw))
      call assigni(nwksz, am(l_nwksz))
      call assigni(iwksz, am(l_iwksz))
      call assigni(nfield, am(l_nfield))

! Call the FMM BEM main program
      call fmmmain(maxa, maxia, am, ia,&
         am(l_n), am(l_x), am(l_y), am(l_node),&
         am(l_dnorm), am(l_bc), am(l_a), am(l_b),&
         am(l_xmax), am(l_xmin), am(l_ymax), am(l_ymin),&
         am(l_ielem), am(l_itree), am(l_level), am(l_loct),&
         am(l_numt), am(l_ifath), am(l_lowlev), am(l_maxl),&
         am(l_levmx), am(l_nexp), am(l_ntylr), am(l_tolerance),&
         am(l_ncellmx), am(l_nleafmx), am(l_mxl), am(l_u),&
         am(l_ax), am(l_nfield), am(l_xfield), am(l_f),&
         am(l_sb), am(l_sx), am(l_igwk), am(l_rgwk),&
         am(l_ligw), am(l_lrgw), am(l_nwksz), am(l_iwksz),&
         am(l_rwork), am(l_iwork))

! Estimate the total CPU time
      call CPU_Time(time)
      write(io_out,*)
      write(*,*)
      write(io_out,*) ' Total CPU time used =', time-time0, '(sec)'
      write(*,*) ' Total CPU time used =', time-time0, '(sec)'
      stop
   end if

end program fmm

!----------------------------------------------------------------------------
! Definition of Variables:
!
! maxa = maximum size of the array am
! maxia = maximum number of variables allowed
! am = a large array storing the variables for the SLATEC GMRES solver
! ia = an array storing the locations of the variables in the array am
!
! n = number of elements (= number of nodes)
! x = coordinates of the nodes
! y = coordinates of the end points of the elements
! node = element connectivity
! dnorm = normal at each node
! bc = BC type and value
!
! a = multipole expansion moments
! b = local expansion coefficients
! xmax,xmin = maximum and minimum x coordinate
! ymax,ymin = maximum and minimum y coordinate
! ielem = ielem(i) gives the original element number for i-th element in
! the quad-tree structure
! itree = itree(c) gives the cell location of c-th cell within each tree level
! loct = elements included in the c-th cell are listed starting at
! the loct(c)-th place in the array ielem
! numt = numt(c) gives the number of elements included in the c-th cell
! ifath = ifath(c) gives the number of the parent cell of the c-th cell
! level = level l cells start at the level(l)-th cell in the tree
! lowlev = number of the tree levels
!
! maxl = maximum number of elements allowed in a leaf
! levmx = maximum number of levels allowed in the tree structure
! nexp = number of terms in multipole expansion
! ntylr = number of terms in local expansion
! tolerance = GMRES solution convergence tolerance
! ncellmx = maximum number of cells allowed in the tree
! nleafmx = maximum number of leaves allowed in the tree
! mxl = maximum dimension of Krylov subspace (used in GMRES)
!
! u = first stores b vector; then solution vector of system Ax = b
! ax = resulting vector of multiplication Ax
! nfield = number of the field points inside the domain
! xfield = coordinates of the field points inside the domain
! f = values of the potential at the field points inside the domain
!
! The following variables and arrays are used in the SLATEC GMRES solver:
! sb,sx,igwk,rgwk,ligw,lrgw,nwksz,iwksz,rwork,iwork
!
!----------------------------------------------------------------------------
